Shifted forces in molecular dynamics 
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Simulations involving the Lennard- Jones potential usually employ a cut-off at r = 2.5a. This 
paper investigates the possibility of reducing the cut-off. Two different cut-off implementations are 
compared, the standard shifted potential cut-off and the less commonly used shifted forces cut-off. 
The first has correct forces below the cut-off, whereas the shifted forces cut-off modifies Newton's 
equations at all distances. The latter is nevertheless superior; we find that for most purposes realistic 
simulations may be obtained using a shifted-forces cut-off at r = 1.5a, even though the pair force is 
here 30 times larger than at r = 2.5a. 



Molecular dynamics (MD) simulations solve Newton's 
equations of motion by discretizing the time coordinate. 
The time-consuming part of any MD simulation is the 
force calculation. For a system of N particles this is an 
0(N 2 ) process whenever all particles interact. In prac- 
tice the interactions are negligible at long distances, how- 
ever, and for this reason one always introduces a cut-off 
at some interparticle distance r = r c beyond which inter- 
actions are ignored 

The standard Lennard- Jones (LJ) pair potential is 
given by 



u W (r) = 4e [(a/r) 12 - {a/rf 



(1) 



Usually, a cut-off at r c = 2.5a is employed; at this point 
the potential is merely 1.6% of its value at the minimum 
(— e). Although a cut-off makes the force calculation an 
almost O(N) process, this calculation remains the most 
demanding in terms of computer time. 

The present paper investigates the possibility of reduc- 
ing the LJ cut-off below 2.5a without compromising accu- 
racy to any significant extent. Before presenting evidence 
that this is possible, it is important to recall that quanti- 
ties depending explicitly on the free energy are generally 
quite sensitive to how large is the cut-off. Examples in- 
clude the location of the critical point @, the surface 
tension 0, H[ , and the solid-liquid coexistence line [3, [f| . 
For such quantities even a cut-off at 2.5a gives inaccu- 
rate results, and in some cases the cut-off must be larger 
than 6(7 to get reliable results [3[. Note, however, that if 
a simulation gives virtually correct particle distribution, 
the thermodynamics of for instance coexisting phases can 
be accurately calculated by application of standard first- 
order perturbation theory 0. 

This note relates to systems for which the standard 
cut-off at 2.5a gives a satisfactory radial distribution 
function. We compared two cut-off implementations at 
varying cut-off's with the "true" LJ system, the latter 
being defined here by the cut-off r c = 4.5a. One cut-off 



is the standard "truncated and shifted potential" (SP for 
shifted potential), for which the radial force is given [l[ 
by = — u' hJ (r) is the LJ radial force) 



, , x J /Lj(r) if r < r c 
fsp(r) = < ., 

if r > r c 



(2) 



This is referred to as a SP cut-off because it corresponds 
to shifting the potential below the cut-off and putting it 
to zero above, which ensures continuity of the potential 
at r c and avoids an infinite force here. 

The "truncated and shifted forces" cut-off (SF for 
shifted forces) [J, 0| has the force go continuously to zero 
at r c , which is obtained by subtracting a constant term: 



fhj(r) - fhj(r c ) if r < r c 
if r > r r 



(3) 
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This corresponds to the following modification of the po- 
tential: Usp(r) = Ulj(V) - (r - r c )u' u (r c ) - u u (r c ) for 
r < r c , usp(r) — for r > r c . Use of a SF cut-off has 
recently become popular in connection with improved 
methods for simulating systems with Coulomb interac- 
tions [8j. 

We simulated the standard single-component LJ liquid 
at the state point that in dimensionless units has density 
p = 0.85 and temperature T = 1.0 [9j. This is a typical 
moderate-pressure liquid state point [l|, [l(| • Other state 
points were also examined - including several state points 
of the fee crystal, at the liquid-gas interface, at the solid- 
liquid interface, and for a supercooled system - leading in 
all cases to the same overall conclusions. For this reason 
we report below results for just one state point of the LJ 
liquid and one of the Kob- Andersen binary LJ (KABLJ) 
liquid [ll| . 2000 LJ particles were simulated using the 
standard central-difference NVT and NVE algorithms 
(Figs. 2, 3 and 4, 6, respectively); 1000 particles of the 
KABLJ liquid were simulated using the NVT algorithm 
(Fig. 5). 

Figure Q] shows the basics of the LJ system. In the 
upper figure the black curve gives the LJ pair potential 
Uhj(f) and the black dashed curve the radial distribu- 
tion function g{r), which has its maximum close to it's 
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FIG. 1: (a) The Lennard- Jones potential (black full curve) 
and the radial distribution function g(r) (black dashed curve) 
for a system at p — 0.85 and T=1.00 in dimensionless units, 
(b) The radial force, fhj{r) = — itLj(r) (black). At r = 1.5a 
the force is 30 times larger than at r = 2.5a. Also shown is 
the shifted force (SF) for a cut-off at 1.5a (red). 
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FIG. 2: Radial distribution function g(r) for the "true" LJ 
system (black) and two cut-off's at r c = 1.5a. The red curve 
gives results for a SF cut-off, the green curve for a SP cut-off. 
The green dashed curve gives results for a SP cut-off with 
smoothing of the force and its derivative at the cut-off [l2| ; 
this however does not improve the SP results. 
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FIG. 3: Integrated numerical difference f Q ' " \g(r) — go(r)\dr 
of the true radial distribution function, go(r), and g(r) for 
various cut-off distances r c . The red curve gives results for 
the SF cut-off, the green for the SP cut-off. Smoothing a SP 
cut-off [l^] does not improve its accuracy (data not shown). 
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FIG. 4: Energy drift as a function of time for long simulations 
(10 8 time steps with h = 0.005) with a cut-off at 1.5a. The 
red curve gives results for the SF cut-off, the green for the SP 
cut-off, and the green dashed curve for a smoothed SP cut- 
off. Smoothing a SP cut-off stabilizes the algorithm, but the 
fluctuations are still somewhat larger than for a SF cut-off. 
The inset shows the initial part of the simulation. 



minimum. In the lower figure the black curve shows the 
LJ pair force The red curve gives /sf(t) when a 

cut-off at 1.5a is introduced; note that the shifted force 
differs significantly from the true force. 

Figure [5] shows the true pair-distribution function 
(black) and the simulated g(r) for three r c — 1.5a cut- 
off's: SF (red), SP (green), and a smoothed SP cut-off 
ensuring the force and its first derivative go continuously 
to zero at the cut-off [HJ (green dashed curve). The 
curves deviate little, except near the cut-off where the 
smallest errors are found for a SF cut-off (inset). 

In order to systematically compare the SP and SF cut- 
off's we studied the LJ liquid for a range of cut-off's. Fig- 
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ure [3] quantifies the difference between the computed g(r) 

and the true, go(r), by evaluating J„ ' 5<T \g(r) ~ go(r)\dr. 
SF is red, SP is green. SF works better than SP for all 
values of r c above the "WCA" cut-off at the potential 
energy minimum H where SF=SP (r c = 2 1 / 6 a = 1.12a). 
Smoothing a SP cut-off has only marginal effect com- 
pared to not smoothing it (results not shown). Applying 
first-order pertubation theory with the g(r) obtained in 
a simulation with SF cut-off at r c — 1.5a leads to a pres- 
sure that deviates only 1% from the correct value. 

Figure U studies energy drift in long NVE simulations 
for r c = 1.5<7. The SF cut-off (red) exhibits no energy 
drift, whereas SP (green) does. FigureHJalso gives results 
when the force of a SP cut-off is smoothed |12j (green 
dashed curve). This leads to much better energy con- 
servation [l|, but the energy fluctuations are somewhat 
larger than for a SF cut-off. The simulations indicate the 
existence of a hidden invariance in the central-difference 
algorithm for a continuous force field deriving from a 
"shadow Hamiltonian" ,13;]. 

Not only static quantities, but also the dynamics are 
affected little by replacing a 2.5<7 SP cut-off with a 1.5a 
SF cut-off. This is demonstrated in Fig. [5j which com- 
pares these two cut-off's for simulations of the incoher- 
ent intermediate scattering function of the supercooled 
KABLJ liquid jll]. For reference a WCA cut-off simula- 
tion is included (blue dashed curve), which was recently 
shown to be inaccurate despite the fact that the WCA 
radial distribution function is reasonably good for this 
system [l4j]- A SP cut-off at r c — 1.5aAA gives too slow 
dynamics (purple dotted curve). Within the numerical 
uncertainties incoherent scattering functions are identi- 
cal for the "true" system, a SP cut-off at r c = 2.5a aa, 
and a SF cut-off at r c = 1.5a aa- Similar results were 
found for the single-component LJ liquid's dynamics. We 
conclude that a SF cut-off at r c — 1.5a generally works 
well for both statics and dynamics of LJ systems. 

Why does a cut-off, for which forces are modified at all 
distances (SF), work better than when forces are correct 
below the cut-off (SP)? A SF cut-off modifies the pair 
force by adding a constant force for all distances below 
r c ; at the same time SF ensures that the pair force goes 
continuously to zero at r = r c . Apparently, ensuring 
continuity of the force - and thereby that u"(r) does not 
spike artificially at the cut-off - is more important than 
maintaining the correct pair force below the cut-off. How 
large is the change induced by the added constant force 
of the SF cut-off? Figure [5] shows the x-component of 
the force on a typical particle as a function of time (r c = 
1.5a). The black curve gives the true force, the red curve 
the SF force, and the blue curve the SF correction term. 
Although the true and SF individual pair forces differ 
significantly (Fig. [TJ, the difference between true and SF 
total forces is small and stochastic (~3%). This reflects 
an almost cancellation of the correction terms deriving 
from the fact that the nearest neighbors are more or less 
uniformly spread around the particle in question. It was 
recently discussed why adding a linear term (oc r) to a 
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FIG. 5: The AA particle incoherent intermediate scattering 
function for the KABLJ liquid in the highly viscous regime 
(T = 0.45, p = 1.2, q = 7.25). The figure shows results for 
the "true" system (black full curve), a SF cut-off at 1.5<taa 
(red full curve), a SP cut-off at 2.5a aa (green dashed curve), 
a SP cut-off at 1.5&AA (purple dotted curve), and for WCA 
cut-off (blue dashed curve). 

pair potential hardly affects dynamics [i~5| and statistical 
mechanics Hi : For a given particle's interactions with 
its neighbors the linear terms sum to almost a constant 
because, if the particle is moved, some nearest-neighbor 
distances increase and others decrease in such a way that 
their sum is almost constant. 

Figure [6jb) shows details of Fig. Ufa); we here added 
the SP force for the same cut-off (green). Both SP and 
correction terms are discontinuous; they jump whenever 
a particle pair distance passes the cut-off. Altogether, 
Fig. [6] shows that not only does the sum of the constant 
forces on a given particle from its neighbors cancel to a 
high degree, so do the interactions with particles beyond 
the cut-off. The result is that the particle distribution 
is affected little by the long-range attractive forces, a 
fact that lies behind the success of perturbation theory 

I ELI El- 

In summary, when a SF cut-off is used instead of the 
standard SP cut-off, errors are significantly reduced. Our 
simulations suggest that a SF cut-off at 1.5cr may be used 
whenever the standard SP cut-off at 2.5a gives reliable 
results; this applies even though the pair force at r — 1.5a 
is 30 times larger than at r = 2.5c A cut-off at 1.5c is 
large enough to ensure that all interactions within the 
first coordination shell are taken into account (Fig. 
Use of a 1.5<7 SF cut-off instead of a SP cut-off at 2.5a 
leads potentially to a factor of (2.5/1.5) 3 = 4.7 shorter 
simulation time for LJ systems. 



Acknowledgments 

The centre for viscous liquid dynamics "Glass and 
Time" is sponsored by the Danish National Research 



4 



i 1— 



340 




_l L 



J l_ 



100 200 300 400 500 600 700 800 900 1000 
Time steps 




360 



380 400 
Time steps 



420 



440 



Foundation (DNRF). 



FIG. 6: (a) The x-component of the force on a typical particle 
during 1000 time steps. The black curve gives the true force, 
the red curve the force for a SF cut-off at 1.5a. The blue curve 
gives the sum of the x-coordinates of the constant "shift" 
terms of Eq. (|3]). (b) Details after 340 steps. The green curve 
gives the SP force (r c = 1.5a). Only true and SF forces are 
smooth functions of time. 
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